#include<math.h>
#include<fstream.h>

class matrix
{
	int row,col;
	double val[100][100];

	public :

	matrix(int,int,double);
	matrix operator+(const matrix &) const;
	matrix operator-(const matrix &) const;
	matrix operator*(const matrix &) const;
	matrix gaussian(const matrix &) const;
	void identity();
	matrix crout(void) const;
	void initialise(double );
	matrix doolittle(void) const;
	matrix choleski(void) const;
	matrix adjoint(void) const;
	int symmetric(void) const;
	matrix inverse_reduction(void) const;
	matrix inverse_adjoint(void) const;
	matrix cofactor(void) const;
	matrix minor_(int,int) const;
	int dominant(int *) const;
	void generate(const int );
	void const_matrix(const matrix &);
	matrix transpose(void) const;
	double determinant(void) const;
	void gauss_jacobi(const matrix &) const;
	void gauss_seidal(const matrix &) const;
	void jacobi_rotate(void) const;
	double power(void) const;
	double inverse_power(void) const;
	void shifted_power(void) const;
	void shifted_inverse_power(void) const;
	double rootmeansqrt(void) const;
	void makezero(matrix &,int,int);
	void makezero(int,int);
	void gerschgorin(double *,int &) const;
	void givens(void) const;
	void householder(void) const;
	double find_s(int &) const;
	void vector_w(double *,double &,int &) const;
	void vector_v(double *,double *,int &) const;
	void vector_q(double *,double *,double *,int &) const;
	void matrix_a(double *,double *,double &,int &);
	double max(void) const;
	double maxelement(int &,int &) const;
	double dist(const matrix &) const;
	double f(int,const double) const;
	void checktol(double tol);
	void sturm_sequence(ofstream &,double *,const int &) const;
	void bisection(ofstream &,double *,double *,const int &) const;
	friend void sort(const double *,const int &);
	friend istream & operator>>(istream &,matrix &);
	friend ostream & operator<<(ostream &,const matrix &);
	friend matrix operator*(const double &,const matrix &);
	friend void gen_matrix(matrix &,matrix &);
	friend void gen_matrix(matrix &);
	friend void write_file(ofstream &,const matrix &);
	friend void write_file(char *,const matrix &,const matrix &,const matrix &);
	friend void write_file(char *,const matrix &,const matrix &);
	friend void write_file(char *,const matrix &,const matrix &,const matrix &,double);
	friend void write_file(char *,const matrix &,const matrix &,double);
	friend void write_file(char *,const matrix &,double);
	friend void write_file(char *,double);
};

matrix::matrix(int x=3,int y=3,double no=0.0)
{	
	row=x;
	col=y;
	
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			val[i][j]=no;


}

void matrix::initialise(double no=0.0)
{
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			val[i][j]=no;


}

void matrix::identity()
{
	for(int i=0;i<row;i++)
		val[i][i]=1;
}

matrix matrix::operator+(const matrix &b) const
{
	matrix c(row,col);

	for(int i=0;i<b.row;i++)
		for(int j=0;j<b.col;j++)
			c.val[i][j]=val[i][j]+b.val[i][j];
	
	return c;
}

matrix matrix::operator-(const matrix &b) const
{
	matrix c(row,col);

	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			c.val[i][j]=val[i][j]-b.val[i][j];

	return c;
}

matrix matrix::operator*(const matrix &b) const
{
	matrix c(row,b.col);
	
	for(int i=0;i<row;i++)
		for(int j=0;j<b.col;j++)
		{
			c.val[i][j]=0;
			for(int k=0;k<col;k++)
				c.val[i][j]=c.val[i][j]+val[i][k]*b.val[k][j];
		}
	
	return c;
}

void matrix::checktol(double tol=0.0001)
{
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			if(fabs(val[i][j])<=tol)
				val[i][j]=0.0;
}

void gen_matrix(matrix &a,matrix &b)
{
	int ch;

	//cout<<"\n Do you want to generate matrix";
	//cout<<"\n (0:Manual 1:Circular 2:Diagonally Dominant 3:For Power 4:General 5:File) : ";
	//cin>>ch;
	ch = 5;
	
	if(ch!=0)
	{
		a.generate(ch);
		
		if(b.col==1 && ch!=5)
			b.const_matrix(a);
		else
			b.generate(ch);

		cout<<"\n Matrix A...\n";
		cout<<a;

		cout<<"\n Matrix B...\n";
		cout<<b;
	}

	else
	{
		cout<<"\n Enter Matrix A of order ("<<a.row<<"X"<<a.col<<") :\n";
		cin>>a;
	
		cout<<"\n Enter Matrix B of order ("<<b.row<<"X"<<b.col<<") :\n";
		cin>>b;
	}
}
	

void gen_matrix(matrix &a)
{
	int ch;

	//cout<<"\n Do you want to generate matrix";
	//cout<<"\n (0:Manual 1:Circular 2:Diagonally Dominant 3:For Power 4:General 5:File) : ";
	//cin>>ch;
	ch=5;
	
	if(ch!=0)
	{
		a.generate(ch);
		
		cout<<"\n Matrix A...\n";
		cout<<a;
	}

	else
	{
		cout<<"\n Enter Matrix A : \n";
		cin>>a;
	}
}

void write_file(ofstream &fout,const matrix &a)
{
	for(int i=0;i<a.row;i++)
	{
		for(int j=0;j<a.col;j++)
			fout<<' '<<a.val[i][j];
			
		fout<<'\n';
	}
}

void write_file(const char *fname,const matrix &a,const matrix &b,const matrix &c)
{
	ofstream fout(fname);
	
	fout<<"\n Matrix1...\n";
	write_file(fout,a);
	fout<<"\n Matrix2...\n";
	write_file(fout,b);
	fout<<"\n Matrix3...\n";
	write_file(fout,c);

}

void write_file(const char *fname,const matrix &a,const matrix &b)
{
	ofstream fout(fname);
	
	fout<<"\n Input Matrix...\n";
	write_file(fout,a);
	fout<<"\n Output Matrix...\n";
	write_file(fout,b);

}

void write_file(const char*fname,const matrix &a,const matrix &b,const double lambda)
{
	ofstream fout(fname);
	
	fout<<"\n Input Matrix...\n";
	write_file(fout,a);
	fout<<"\n Vector matrix...\n";
	write_file(fout,b);
	fout<<"\n Lambda...\n";
	fout<<lambda;
}

void write_file(const char*fname,const matrix &a,const matrix &b,const matrix &c,const double d)
{
	ofstream fout(fname);
	
	fout<<"\n Input Matrix...\n";
	write_file(fout,a);
	fout<<"\n Lower Triangular Matrix...\n";
	write_file(fout,b);
	fout<<"\n Upper Triangular Matrix...\n";
	write_file(fout,c);
	fout<<"\n Determinant...\n";
	fout<<d;
}

void write_file(const char*fname, const matrix &a, const double d)
{
	ofstream fout(fname);
	
	fout<<"\n Input Matrix...\n";
	write_file(fout,a);
	fout<<"\n Value...\n";
	fout<<d;
}

matrix matrix::transpose(void) const
{
	matrix b(row,col);
	for(int i=0;i<col;i++)
		for(int j=0;j<row;j++)
			b.val[i][j]=val[j][i];
	
	return b;
}

matrix matrix::minor_(int i,int j) const
{
	matrix submat(row-1,col-1);
	int m=0,n=0;

	for(int k=0;k<row;k++)
	{
		if(k!=i)
		{
			for(int l=0;l<col;l++)
				if(l!=j)
					submat.val[m][n++]=this->val[k][l];
		
			n=0;
			m++;
		}
	}
	
	return submat;
}

matrix matrix::cofactor(void) const
{
	
	matrix cofact(row,col),submat(row-1,col-1);
	
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
		{
			submat=this->minor_(i,j);
			
			if(((i+j)%2)==1)
				cofact.val[i][j]=-1*submat.determinant();
			else
				cofact.val[i][j]=submat.determinant();

		}
	
	return cofact;
}


double matrix::determinant(void) const
{
	double det=1.0,m;
	int record[row],temp;
	matrix b(row,col);
	
	b=*this;
	
	/*Initialize row value*/
	for(int i=0;i<row;i++)
		record[i]=i;
	
	/*Process to find Upper Triangular*/
	for(int i=0;i<(row-1);i++)
	{
		for(int j=i+1;j<row;j++)
		{
			if(fabs(b.val[record[j]][i])>fabs(b.val[record[i]][i]))	/*Pivoting process*/
			{
				det=(-1)*det;
				temp=record[i];
				record[i]=record[j];
				record[j]=temp;
			}
			
		}

		if(fabs(b.val[record[i]][i])<0.0001)
		{
				b.val[record[i]][i]=0;
				det=0;
				return det;
		}
		
		for(int j=i+1;j<row;j++)
		{
			m=b.val[record[j]][i]/b.val[record[i]][i];	/*Find ratio*/
			
			b.val[record[j]][i]=0;
			for(int k=i+1;k<col;k++)
			{
				b.val[record[j]][k]=b.val[record[j]][k]-m*b.val[record[i]][k];
				
				if(fabs(b.val[record[j]][k])<=0.0001)
					b.val[record[j]][k]=0;
			}
		}
	}
	
	for(int i=0;i<row;i++)
		det=det*b.val[record[i]][i];
		
	write_file("../out/determinant.txt",*this,det);
	
	return det;
}

matrix matrix::inverse_reduction(void) const
{
	double m;
	matrix invert(row,col),b(row,col),dum(row,col);
	int record[row],temp;
	
	b=*this;
	
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			if(i==j)
				invert.val[i][j]=1;
			else
				invert.val[i][j]=0;

	
	/*Initialize row value*/
	for(int i=0;i<row;i++)
		record[i]=i;
		
	
	/*Process to find Identity Matrix*/
	for(int i=0;i<row;i++)
	{
		if((b.val[record[i]][i]!=1.0)&&fabs(b.val[record[i]][i])>0.0001)
		{
			for(int j=0;j<row;j++)
			{
				
				if(j!=i)
					b.val[record[i]][j]=b.val[record[i]][j]/b.val[record[i]][i];
					
				invert.val[record[i]][j]=invert.val[record[i]][j]/b.val[record[i]][i];
			}

			//invert.val[record[i]][i]=invert.val[record[i]][i]/val[record[i]][i];
			b.val[record[i]][i]=1;
		}
		
		else if(fabs(b.val[record[i]][i])<=0.0001)
		{
			for(int j=i+1;j<row;j++)
			{
				if(fabs(b.val[record[j]][i])>fabs(b.val[record[i]][i]))
				{
					temp=record[i];
					record[i]=record[j];
					record[j]=temp;
				}
			}
		}
		
		b.checktol();
		for(int j=0;j<row;j++)
		{
			
			if(i!=j)
			{
				m=b.val[record[j]][i];	/*Find ratio*/
				
				b.val[record[j]][i]=0;
				
				for(int k=0;k<col;k++)
				{
						if(k!=i)
							b.val[record[j]][k]=b.val[record[j]][k]-m*b.val[record[i]][k];
						invert.val[record[j]][k]=invert.val[record[j]][k]-m*invert.val[record[i]][k];
						
				}

			}
		}

	}

	cout<<"\n Inverse...\n";
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			dum.val[i][j]=invert.val[record[i]][j];
			
	write_file("../out/inverse_reduction.txt",*this,dum);

	return dum;

}

matrix matrix::crout(void) const
{
	matrix l(row,col),u(row,col),lu(row,col);
	matrix l1(row,col),u1(row,col);
	double sum,det=1;
	int record[row];

	for(int i=0;i<row;i++)
		record[i]=i;

	up:for(int i=0;i<row;i++)
	{

		for(int j=0;j<col;j++)
			if(i<j)
			{
				sum=0;
				for(int k=0;k<i;k++)
					sum=sum+l.val[record[i]][k]*u.val[record[k]][j];



				if(fabs(l.val[record[i]][i])<=0.0001)
				{
					if(i==(row-1))
					{
						int temp=record[i];
						record[i]=record[i-1];
						record[i-1]=temp;
						i=0;
						j=0;
						det=-det;
						l.initialise(0);
						u.initialise(0);
						goto up;
					}
					else
					{
						int temp=record[i];
						record[i]=record[i+1];
						record[i+1]=temp;
						i=0;
						j=0;
						l.initialise(0);
						u.initialise(0);
						det=-det;
						goto up;
					}
				}
				
				u.val[record[i]][j]=(1/l.val[record[i]][i])*(val[record[i]][j]-sum);
				
			}
			
			else if(i>=j)
			{
				if(i==j)
					u.val[record[i]][j]=1;
								
				sum=0;
				for(int k=0;k<j;k++)
					sum=sum+l.val[record[i]][k]*u.val[record[k]][j];
				
				l.val[record[i]][j]=val[record[i]][j]-sum;
				
				
			}
	}
	
	
	
	
	for(int i=0;i<row;i++)
	{
		for(int j=0;j<col;j++)
		{
			l1.val[i][j]=l.val[record[i]][j];
			
			if(i==j)
				det=det*l.val[record[i]][j];
			
		}
		cout<<endl;
	}
			

	for(int i=0;i<row;i++)
	{
		for(int j=0;j<col;j++)
			u1.val[i][j]=u.val[record[i]][j];

		cout<<endl;
	}
	
	cout<<"\n L Matrix...\n";
	cout<<l1;
	
	cout<<"\n U Matrix...\n";
	cout<<u1;
	
	cout<<"\n Determinant(A) = "<<det;
	
	lu=l1*u1;
	
	write_file("../out/crout.txt",*this,l1,u1,det);
	return(lu);
}

matrix matrix::doolittle(void) const
{
	matrix l(row,col),u(row,col),lu(row,col);
	matrix l1(row,col),u1(row,col);
	double sum,det=1;
	int record[row];
	
	for(int i=0;i<row;i++)
		record[i]=i;
	
	up:for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			if(i<=j)
			{	
				if(i==j)
					l.val[record[i]][j]=1;

				sum=0;
				for(int k=0;k<j;k++)
					sum=sum+l.val[record[i]][k]*u.val[record[k]][j];

				u.val[record[i]][j]=val[record[i]][j]-sum;
				
			}
			else if(i>j)
			{
				sum=0;
				for(int k=0;k<i;k++)
					sum=sum+l.val[record[i]][k]*u.val[record[k]][j];
				
				if(fabs(u.val[record[j]][j])<=0.0001)
				{
					if(j==(row-1))
					{
						int temp=record[j];
						record[j]=record[j-1];
						record[j-1]=temp;
						i=0;
						j=0;
						det=-det;
						l.initialise(0);
						u.initialise(0);
						goto up;
					}
					else
					{
						int temp=record[j];
						record[j]=record[j+1];
						record[j+1]=temp;
						i=0;
						j=0;
						l.initialise(0);
						u.initialise(0);
						det=-det;
						goto up;
					}
				}
				l.val[record[i]][j]=(1/u.val[record[j]][j])*(val[record[i]][j]-sum);
			}
		
	for(int i=0;i<row;i++)
	{
		for(int j=0;j<col;j++)
			l1.val[i][j]=l.val[record[i]][j];
	}
			

	for(int i=0;i<row;i++)
	{
		for(int j=0;j<col;j++)
		{
			u1.val[i][j]=u.val[record[i]][j];
		
			if(i==j)
				det=det*u.val[record[i]][j];
		}
	}
	
	cout<<"\n L Matrix...\n";
	cout<<l1;
	
	cout<<"\n U Matrix...\n";
	cout<<u1;

	cout<<"\n Determinant(A) = "<<det;
	
	lu=l1*u1;
	write_file("../out/doolittle.txt",*this,l1,u1,det);
	return(lu);
}

matrix matrix::choleski(void) const
{
	matrix l(row,col),u(row,col),lu(row,col);
	double sum,det=1;
	
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			if(i>=j)
			{
				if(i==j)
				{
					sum=0;
					for(int k=0;k<i;k++)
						sum=sum+l.val[i][k]*l.val[i][k];
					
					l.val[i][i]=sqrtl(fabs(val[i][i]-sum));
					det=det*l.val[i][j];
				}
				else
				{
					sum=0;
					for(int k=0;k<i;k++)
						sum=sum+l.val[i][k]*l.val[j][k];

					l.val[i][j]=(1/l.val[j][j])*(val[i][j]-sum);
				
				}
			}
				
	
	cout<<"\n L Matrix...\n";
	cout<<l;
	
	u=l.transpose();
	cout<<"\n Transpose of L Matrix...\n";
	cout<<u;
	
	cout<<det*det;
	
	cout<<"\n Determinant(A) = "<<det;

	lu=l*u;

	write_file("../out/choleski.txt",*this,l,u,det);
	return(lu);
}


matrix matrix::inverse_adjoint(void) const
{
	double det;
	matrix invert(row,col),adj(row,col);

	det=this->determinant();

	adj=this->adjoint();

	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			invert.val[i][j]=(1./det)*adj.val[i][j];

	write_file("../out/inverse_adj.txt",*this,invert);
	return invert;
}


	
istream & operator>>(istream &din, matrix &b)
{
	for(int i=0;i<b.row;i++)
		for(int j=0;j<b.col;j++)
			din>>b.val[i][j];
	
	return din;
}

ostream & operator<<(ostream &dout, const matrix &b)
{
	for(int i=0;i<b.row;i++)
	{
		for(int j=0;j<b.col;j++)
		{
			dout.width(12);
			dout<<b.val[i][j];
		}
	
		dout<<'\n';
	}
	
	return dout;
}

void matrix::generate(const int x)
{
	if(x==1)
	{
		for(int i=0;i<row;i++)
			for(int j=0;j<col;j++)
				if(abs(i-j)==1)
					val[i][j]=-10;
				else 
					val[i][j]=0;

		val[0][col-1]=val[row-1][0]=-10;
	}
	
	else if(x==2)
	{
		double sum;

		for(int i=0;i<row;i++)
	        {
                	sum=0.0;
	                for(int j=0;j<col;j++)
        	        {
                	        if(i!=j)
                        	{
                                	val[i][j]=(double) (10.0*rand()/(RAND_MAX+1.0));
	                                sum=sum+fabs(val[i][j]);
        	                }
                	}

	                if(sum>=100)
        	                val[i][i]=sum+1+(double) (10.0*rand()/(RAND_MAX+1.0));
                	else
	        		val[i][i]=100+(double) (10.0*rand()/(RAND_MAX+1.0));

         	}
	}

	else if(x==3)
		for(int i=0;i<row;i++)
			for(int j=0;j<col;j++)
				if(i==j)
					val[i][j]=2;
				else
					val[i][j]=1;
	else if(x==4)
		for(int i=0;i<row;i++)
			for(int j=0;j<col;j++)
				val[i][j]=1+(double) (10.0*rand()/(RAND_MAX+1.0));
	else if(x==5)
	{

		char name[15];
		cout<<"\n Enter the file name : ";
		cin>>name;
		
		char path[50]="../inp/";
		ifstream fin(strcat(path,name));
	 	
		if(fin.fail())
		{
			cout<<"\n Could not open the file...";
			return;
		}

		
		if(strcmp(name,"const.txt")!=0)
		{
			fin>>row;
			col=row;
		}
		else
			fin>>row;
			
		for(int i=0;i<row;i++)
			for(int j=0;j<col;j++)
				if(fin)
					fin>>val[i][j];	
		fin.close();
	}
}

void matrix::const_matrix(const matrix &b)
{
	double sum;

	for(int i=0;i<row;i++)
	{
		sum=0.0;
		for(int j=0;j<row;j++)
			sum=sum+b.val[i][j];
		
		val[i][0]=sum;
	}
}

matrix matrix::gaussian(const matrix &value) const
{
	double m;
	matrix x(row,1);
	int record[row],temp;
	matrix b(row,col),c(row,1);
	b=*this;
	c=value;

	/*Initialize row value*/
	for(int i=0;i<row;i++)
		record[i]=i;

	/*Process to find Upper Triangular*/
	for(int i=0;i<(row-1);i++)
	{
		for(int j=i+1;j<row;j++)
		{
			if(fabs(b.val[record[j]][i])>fabs(b.val[record[i]][i]))	/*Pivoting process*/
			{
				temp=record[i];
				record[i]=record[j];
				record[j]=temp;
			}
		}

		for(int j=i+1;j<row;j++)
		{
			m=b.val[record[j]][i]/b.val[record[i]][i];	/*Form multiplier*/

			b.val[record[j]][i]=0;

			for(int k=i+1;k<row;k++)
				b.val[record[j]][k]=b.val[record[j]][k]-m*b.val[record[i]][k];/*End UT*/

			c.val[record[j]][0]=c.val[record[j]][0]-m*c.val[record[i]][0];
		}
	}

	/*Process of the back substitution*/
	double sum;

	x.val[row-1][0]=c.val[record[row-1]][0]/b.val[record[row-1]][row-1];

	for(int i=row-2;i>=0;i--)
	{
		sum=0;
		for(int j=i+1;j<col;j++)
			sum=sum+b.val[record[i]][j]*x.val[j][0];


		x.val[i][0]=(c.val[record[i]][0]-sum)/b.val[record[i]][i];
	}

	write_file("../out/guassian.txt",b,c,x);

	return x;
}

int matrix::symmetric(void) const
{

	for(int i=0;i<row;i++)
		for(int j=0;j<i;j++)
			if(val[i][j]!=val[j][i])
				return 0;

	return 1;
}


matrix matrix::adjoint(void) const
{

	matrix adj(row,col);

	adj=this->cofactor();

	return adj.transpose();
}

matrix operator*(const double &n,const matrix &b)
{
	matrix c(b.row,b.col);

	for(int i=0;i<b.row;i++)
		for(int j=0;j<b.col;j++)
			c.val[i][j]=n*b.val[i][j];

	
	return c;
}


int matrix::dominant(int *rec) const
{
	int flag;
	double sum;
	for(int i=0;i<row;i++)
	{
		flag=0;
		sum=0;
		for(int j=0;j<col;j++)
			 sum=sum+fabs(val[i][j]);

		for(int j=0;j<col;j++)
			if(2*fabs(val[i][j])>sum)
			{
				rec[j]=i;
				flag=1;

				for(int k=0;k<j;k++)
					if(rec[k]==i)
						flag=0;

				break;
			}
	}

        for(int i=0;i<row;i++)
        	if(rec[i]==-1)
                	flag=0;

	return flag;
}



void matrix::gauss_jacobi(const matrix & c) const
{
	matrix x(row,1),xnew(row,1);
	int record[row],count,iteration,d;
	double sum,error,norm;

	cout<<"\n Enter the error & no. of iteration : ";
	cin>>error>>iteration;

	/*Initialize row value*/
	for(int i=0;i<row;i++)
		record[i]=-1;

	d=this->dominant(record);

	cout<<"\n Coefficient Matrix  A ...\n";
        cout<<*this;
	cout<<"\n Constant Matrix B ...\n";
        cout<<c;

        if(d==0)
	{
		cout<<"\n Matrix is not diagonally dominant\n";
		return;
	}
        else
        	cout<<"\n Matrix is diagonally dominant\n";

	for(int i=0;i<row;i++)
	{
		x.val[record[i]][0]=c.val[record[i]][0]/val[record[i]][i];
		
		if(fabs(x.val[record[i]][0])<=error)
			x.val[record[i]][0]=0;
	}
	
	for(count=1;count<=iteration;count++)
	{
		for(int i=0;i<row;i++)
		{
			sum=c.val[record[i]][0];

			for(int j=0;j<col;j++)
				if(i!=j)
					sum-=val[record[i]][j]*x.val[record[j]][0];

			xnew.val[record[i]][0]=sum/val[record[i]][i];
		
			if(fabs(x.val[record[i]][0])<=error)
				xnew.val[record[i]][0]=0;
		}


		norm=x.dist(xnew);

                for(int j=0;j<col;j++)
			x.val[record[j]][0]=xnew.val[record[j]][0];

		if(norm<=error)
			break;
	}


	cout<<"\n Total no. of iteration required = "<<count;
	cout<<"\n Variable Matrix X...\n";
	cout<<xnew;
	
	write_file("../out/gauss_jacobi.txt",*this,c,xnew);

}

void matrix::gauss_seidal(const matrix & c) const
{
	matrix x(row,1),xnew(row,1);
	int record[row],count,iteration,d;
	double sum,error,norm;

	cout<<"\n Enter the error & no. of iteration : ";
	cin>>error>>iteration;

	/*Initialize row value*/
	for(int i=0;i<row;i++)
		record[i]=-1;

	d=dominant(record);

	cout<<"\n Coefficient Matrix  A ...\n";
        cout<<*this;
	cout<<"\n Constant Matrix B ...\n";
        cout<<c;

	if(d==0)
	{
		cout<<"\n Matrix is not diagonally dominant\n";
		return;
	}
        else
        	cout<<"\n Matrix is diagonally dominant\n";



	for(int i=0;i<row;i++)
	{
		double temp;
		temp=c.val[record[i]][0]/val[record[i]][i];

		if(fabs(temp)<=error)
			temp=0;

		x.val[record[i]][0]=temp;

		xnew.val[record[i]][0]=temp;
	}


	for(count=1;count<=iteration;count++)
	{
		for(int i=0;i<row;i++)
		{
			sum=c.val[record[i]][0];

			for(int j=0;j<col;j++)
				if(i!=j)
					sum-=val[record[i]][j]*xnew.val[record[j]][0];

			xnew.val[record[i]][0]=sum/val[record[i]][i];

			if(fabs(xnew.val[record[i]][0])<=error)
				xnew.val[record[i]][0]=0;

		}


		norm=x.dist(xnew);

                for(int j=0;j<col;j++)
			x.val[record[j]][0]=xnew.val[record[j]][0];

		if(norm<=error)
			break;
	}

	cout<<"\n Total no. of iteration required = "<<count;
	cout<<"\n Variable Matrix X...\n";
	cout<<xnew;

	write_file("../out/gauss_seidal.txt",*this,c,xnew);
}

double matrix::max(void) const
{

	double max=fabs(val[0][0]),c=val[0][0];
	for(int i=0;i<row;i++)
			if(fabs(val[i][0])>max)
			{
				max=fabs(val[i][0]);
				c=val[i][0];
			}

	return c;
}

double matrix::maxelement(int &m,int &n) const
{

	double max=fabs(val[0][1]);
	m=0;
	n=1;
	for(int i=0;i<row-1;i++)
		for(int j=i+1;j<row;j++)
			if(fabs(val[i][j])>max)
			{
				max=fabs(val[i][j]);
				m=i;
				n=j;

			}

	return max;

}
double matrix::dist(const matrix &b) const
{
	double d=0.0;

	for(int i=0;i<row;i++)
		d=d+(val[i][0]-b.val[i][0])*(val[i][0]-b.val[i][0]);

	return sqrtl(d);
}

double matrix::power(void) const
{
	double lambda=0,c,error,diff,d;
	int iteration;
	matrix b(row,col),x(row,1,1.0),y(row,1);

	b=*this;

	cout<<"\n Enter error & total iteration : ";
	cin>>error>>iteration;

	cout<<"\n x = \n";
	cout<<x;
	for(int i=1;i<=iteration;i++)
	{
		y=b*x;

		c=y.max();
		y=(1/c)*y;

		y.checktol();

		diff=fabs(lambda-c);
		lambda=c;
		d=y.dist(x);

		x=y;

		if(((diff>d)?diff:d)<error)
		{
			iteration=i;
			break;
		}
	}

	cout<<"\n Total iteration required : "<<iteration;
	cout<<"\n Eigenpair of A is ...";
	cout<<"\n lambda = "<<lambda;

	cout<<"\n vector x = \n";
	cout<<y;

	write_file("../out/power.txt",*this,y,lambda);
	return lambda;
}

double matrix::inverse_power(void) const
{
	double lambda=0,c,error,diff,d;
	int iteration;
	matrix b(row,col),x(row,1,1.0),y(row,1);

	b=*this;

	cout<<"\n Enter error & total iteration : ";
	cin>>error>>iteration;
	
	cout<<"\n x = \n";
	cout<<x;
	for(int i=1;i<=iteration;i++)
	{
		y=b.gaussian(x);

		c=y.max();
		y=(1/c)*y;
		
		y.checktol();
		diff=fabs(lambda-c);
		lambda=c;
		d=y.dist(x);
		
		x=y;
		
		if(((diff>d)?diff:d)<error) 
		{
			iteration=i;
			break;
		}
	}
	
	cout<<"\n Total iteration required : "<<iteration;
	cout<<"\n Eigenpair of A is ...";
	cout<<"\n lambda = "<<(1/lambda);
	
	cout<<"\n Eigenpair of inverse(A) is ...";
	cout<<"\n lambda = "<<lambda;
	cout<<"\n vector x = \n";
	cout<<x;

	write_file("../out/inverse_power.txt",*this,y,lambda);
	return lambda;
}

void matrix::shifted_power() const
{
	double lambda=0,c,error,diff,d,s;
	int iteration;
	matrix b(row,col),x(row,1,1.0),y(row,1),I(row,col);

	b=*this;

	cout<<"\n Enter error & total iteration : ";
	cin>>error>>iteration;

	cout<<"\n Enter shifting element : ";
	cin>>s;
	
	for(int i=0;i<row;i++)
		I.val[i][i]=1;
		
	b=b-(s*I);
	for(int i=1;i<=iteration;i++)
	{
		y=b*x;
		
		c=y.max();
		y=(1/c)*y;
		
		y.checktol();
		diff=fabs(lambda-c);
		lambda=c;
		d=y.dist(x);

		x=y;

		if(((diff>d)?diff:d)<error)
		{
			iteration=i;
			break;
		}
	}

	cout<<"\n Total iteration required : "<<iteration;
	cout<<"\n Eigenvalue of A is ...";
	cout<<"\n lambda = "<<lambda+s;
	cout<<"\n Eigenpair of (A-sI) is ...";
	cout<<"\n lambda = "<<lambda;
	cout<<"\n vector x = \n";
	cout<<x;
	
	write_file("../out/shift.txt",*this,y,lambda);
}

void matrix::shifted_inverse_power() const
{
	double lambda=0,c,error,diff,d,s;
	int iteration;
	matrix b(row,col),x(row,1,1.0),y(row,1),I(row,col);

	b=*this;

	cout<<"\n Enter error & total iteration : ";
	cin>>error>>iteration;

	cout<<"\n Enter shifting element : ";
	cin>>s;

	for(int i=0;i<row;i++)
		I.val[i][i]=1;

	b=b-(s*I);
	
	for(int i=1;i<=iteration;i++)
	{
		y=b.gaussian(x);
		
		c=y.max();
		y=(1/c)*y;
		
		y.checktol();
		diff=fabs(lambda-c);
		lambda=c;
		d=y.dist(x);

		x=y;
		
		if(((diff>d)?diff:d)<error)
		{
			iteration=i;
			break;
		}
	}
	
	
	cout<<"\n Total iteration required : "<<iteration;
	cout<<"\n Eigevalue of A is ...";
	cout<<"\n lambda = "<<((1/(lambda))+s);
	cout<<"\n Eigenpair of inverse(A-sI) is ...";
	cout<<"\n lambda = "<<lambda;
	cout<<"\n vector x = \n";
	cout<<x;
	
	write_file("../out/shift_inverse.txt",*this,y,lambda);
}


double matrix::rootmeansqrt(void) const
{
	double sum=0.0;
	for(int i=0;i<row;i++)
		sum=sum+val[i][i]*val[i][i];

	return sqrt(sum/row);
}

void matrix::makezero(matrix &v,int p,int q)
{
	double angle,t,c,s,xp[col],xq[col],yp[col],yq[col];
	int sign;

	angle=(val[q][q]-val[p][p])/(2*val[p][q]);

	(angle>=0)?sign=1:sign=-1;

	t=sign/(fabs(angle)+sqrt(angle*angle+1));
	c=1/sqrt(t*t+1);
	s=c*t;
	
	cout<<"\n c = "<<c;
	cout<<"\n s = "<<s;
	
	for(int i=0;i<row;i++)
	{
		xp[i]=val[i][p];
		xq[i]=val[i][q];
	}

	val[p][q]=0;
	val[q][p]=0;
	val[p][p]=c*c*xp[p]+s*s*xq[q]-2*c*s*xq[p];
	val[q][q]=s*s*xp[p]+c*c*xq[q]+2*c*s*xq[p];

	for(int i=0;i<row;i++)
	{
		if((i!=p) && (i!=q))
		{
			val[i][p]=c*xp[i]-s*xq[i];
			val[p][i]=val[i][p];
			
			val[i][q]=c*xq[i]+s*xp[i];
			val[q][i]=val[i][q];
			
		}

		yp[i]=v.val[i][p];
		yq[i]=v.val[i][q];
	}

	for(int i=0;i<row;i++)
	{
		v.val[i][p]=c*yp[i]-s*yq[i];
		v.val[i][q]=s*yp[i]+c*yq[i];
	}


}


void matrix::jacobi_rotate(void) const
{
	double rms,err,max;
	matrix v(row,col),b(row,col);
	int maxs,counts,flag=1,m,n;

	if(this->symmetric()!=1)
	{
		cout<<"\n Matrix is not symmetric";
		return ;
	}

	b=*this;

	for(int i=0;i<row;i++)
		v.val[i][i]=1;

	cout<<"\n Enter tolerance & maximum iterations : ";
	cin>>err>>maxs;

	while((flag==1) && (maxs>=counts))
	{
		rms=b.rootmeansqrt();
		counts=counts+1;
		flag=0;

		max=b.maxelement(m,n);
		cout<<"\n Iteration   : "<<counts;
		cout<<"\n Max Element : "<<max;

		if(max/rms > err)
		{
			b.makezero(v,m,n);
			b.checktol(err);
			flag=1;
		}
	}

	cout<<"\n Total iteration required = "<<counts;
	cout<<"\n Diagonal matrix is...\n";
	cout<<b;

	cout<<"\n The eigen vectors of matrix A are...\n";
	cout<<v;

	cout<<"\n The eigen values of matrix A are...\n";
	for(int i=0;i<row;i++)
		cout<<"\n "<<b.val[i][i];
		
	write_file("../out/jacobi_rotate.txt",*this,v,b);

}

void matrix::gerschgorin(double *lambda,int &c) const
{
	double pk;
	int flag;

	for(int i=0;i<row;i++)
	{
		pk=0.0;
		for(int j=0;j<row;j++)
			if(j!=i)
				pk=pk+fabs(val[i][j]);

		for(double j=-pk+val[i][i];j<=pk+val[i][i];j++)
		{
			flag=1;
			for(int k=0;k<=c;k++)
				if(fabs(lambda[k]-j)<=0.0001)
				{
					flag=0;
					break;
				}

			if(flag)
				lambda[++c]=j;
		}

	}
}

double matrix::f(int n,const double lambda) const
{
	matrix temp(n,n),I(n,n);
	double det;

	for(int i=0;i<n;i++)
		I.val[i][i]=1;

	temp=lambda*I;
	temp=temp-*this;

	det=temp.determinant();

	return det;
}

void matrix::makezero(int p,int q)
{
	double t,c,s,xp[col],xq[col];
	matrix rot(row,col),rott(row,col);

	if(fabs(val[p][p+1])>0.0001)
        	t=val[p][q]/val[p][p+1];
	else
		t=3.14/2;
		
        c=1/sqrt(t*t+1);
	s=c*t;

	cout<<"\n tan = "<<t;
	cout<<"\n c = "<<c;
	cout<<"\n s = "<<s;

	for(int i=0;i<row;i++)
		rot.val[i][i]=1;

	p=p+1;
	
	cout<<"\n Rotation at ("<<p<<", "<<q<<")\n";
	rot.val[p][p]=c;
	rot.val[p][q]=-s;
	rot.val[q][p]=s;
	rot.val[q][q]=c;

	rott=rot.transpose();
	*this=rott*(*this);
	*this=(*this)*rot;

}


void sort(double *lambda, const int &c)
{
	
	for(int i=0;i<=c;i++)
		for(int j=i+1;j<=c;j++)
			if(lambda[i]>lambda[j])
			{
				double temp=lambda[i];
				lambda[i]=lambda[j];
				lambda[j]=temp;
			}

}

void matrix::bisection(ofstream &fout,double *lambda,double *v,const int &c) const
{
	double f1,f2,x0,x1,x2;
	int r=row;
	
	for(int i=0;i<=c;i++)
	{
		r=row;
		x0=0;
		x1=0;

		if(v[i]==-2)
		{
			cout<<"\n Eigenvalue = "<<lambda[i];
			fout<<"\n Eigenvalue = "<<lambda[i];
		}

		else if((i+1)>c)
			continue;
                else if(fabs(v[i]-v[i+1])!=0 && (v[i+1]!=-2))
		{
			double f1,f2;
			
			up:f1=f(r,lambda[i]);
			f2=f(r,lambda[i+1]);
			
			if(f1<0 && f2>=0)
			{
				x0=lambda[i];
				x1=lambda[i+1];
			}
			else if(f1>=0 && f2<0)
			{
				x0=lambda[i+1];
				x1=lambda[i];
			}
			
			else
			{
				r--;
				if(r<0)
					continue;
				goto up;
			}
				
			while(fabs(x1-x0)>0.0001)
			{
                		x2=(x0+x1)/2;

				if(f(row,x2)<0)
        	                	x0=x2;
                		else
					x1=x2;
			}

			cout<<"\n Eigenvalue = "<<x2;
			fout<<"\n Eigenvalue = "<<x2;

		}
	}
		
}

void matrix::sturm_sequence(ofstream &fout,double *lambda,const int &c) const
{
	
	double v[c+1],fval,x0,x1,f1,f2;
	int no=c,t;
	
	sort(lambda,c);
	
	cout<<"\n Sturm sequences as ...\n";
	
	cout<<"\n\tlambda";
	for(int i=0;i<=row;i++)
		cout<<"\tf("<<i<<")";
	cout<<"\tv(lambda)";
	

	for(int i=0;i<=no;i++)
	{
		t=1;
		v[i]=0;
		cout<<"\n\t"<<lambda[i];
		
		for(int j=0;j<=row;j++)
		{
			if(j==0)
				t=1;
			
			else
			{
				fval=f(j,lambda[i]);

				if(fabs(fval)<=0.0001)
					fval=0;
				
				if(fabs(fval)<=0.0001 && j==row)
				{
					v[i]=-2;
					t=1;
				}
				
				else if(fval>=0 && t==-1)
				{
					v[i]=v[i]+1;
					t=1;
				}
				
				else if(fval<0 && t==1)
				{
					v[i]=v[i]+1;
					t=-1;
				}
				
				else if(fval<0)
					t=-1;
				
				else if(fval>=0)
					t=1;
			}
			
			if(t==1)
				cout<<"\t+";
			else
				cout<<"\t-";
		}
		
		if(v[i]==-2)
			cout<<"\tx";
		else
			cout<<"\t"<<v[i];
	}

	
	this->bisection(fout,lambda,v,c);
	
}

void matrix::householder(void) const
{
	matrix b(row,col);
	double w[row],v1[row],lambda[100],s,q[100];
	int counts=0,c=-1,v[100];
	ofstream fout("../out/householder.txt");

	if(this->symmetric()!=1)
	{
		cout<<"\n Matrix is not symmetric";
		return ;
	}

	b=*this;
	
	for(int k=0;k<row-2;k++)
	{
		cout<<"\n Iteration = "<<k+1;
		s=b.find_s(k);
		b.vector_w(w,s,k);
		b.vector_v(v1,w,k);
		b.vector_q(q,v1,w,k);
		b.matrix_a(q,w,s,k);
	}		
	
	cout<<"\n Tridiagonal Matrix...\n";
	cout<<b;

	this->gerschgorin(lambda,c);

	b.sturm_sequence(fout,lambda,c);
	
	fout<<"\n Input Matrix...\n";
	write_file(fout,*this);
	fout<<"\n Triangular Matrix...\n";
	write_file(fout,b);
	
}

double matrix::find_s(int &k) const
{
	double s=0;
	for(int i=k+1;i<row;i++)
		s=s+val[i][k]*val[i][k];

	s=sqrtl(s);

	if(val[k+1][k]<0)
		s=-s;

	//cout<<"\n S = "<<s;
	return s;
}

void matrix::vector_w(double *w,double &s,int &k) const
{
	double r;


	r=sqrtl(2*val[k+1][k]*s+2*s*s);
	
	for(int i=0;i<=k;i++)
		w[i]=0;

	w[k+1]=(val[k+1][k]+s)/r;

	for(int i=k+2;i<row;i++)
		w[i]=val[i][k]/r;

	//cout<<"\n R = "<<r;
	//cout<<"\n WT = ";
	//for(int i=0;i<row;i++)
	//	cout<<' '<<w[i];


}
		
void matrix::vector_v(double *v,double *w,int &k) const
{
	double sum;

	for(int i=0;i<=k;i++)
		v[i]=0;

	for(int i=k+1;i<row;i++)
	{
		sum=0;

		for(int j=k+1;j<row;j++)
			sum=sum+val[i][j]*w[j];

		v[i]=sum;
	}

	//cout<<"\n VT = ";
	//for(int i=0;i<row;i++)
	//	cout<<' '<<v[i];

}

void matrix::vector_q(double *q,double *v,double *w,int &k) const
{
	double c=0;
	for(int i=k+1;i<row;i++)
		c=c+v[i]*w[i];
	
	for(int i=0;i<=k;i++)
		q[i]=0;

	for(int i=k+1;i<row;i++)
		q[i]=v[i]-c*w[i];

	//cout<<"\n c = "<<c;

	//cout<<"\n QT = ";
	//for(int i=0;i<row;i++)
	//	cout<<' '<<q[i];
}

void matrix::matrix_a(double *q,double *w,double &s,int &k)
{
	for(int i=k+2;i<row;i++)
		val[i][k]=val[k][i]=0;

	val[k+1][k]=val[k][k+1]=-s;

	for(int i=k;i<row;i++)
		val[i][i]=val[i][i]-4*q[i]*w[i];

	for(int i=k+1;i<row;i++)
		for(int j=i+1;j<row;j++)
			val[i][j]=val[j][i]=val[i][j]-2*w[i]*q[j]-2*q[i]*w[j];

	checktol();
	cout<<"\n Matrix A = \n";
	cout<<*this;
			
}

		

void matrix::givens(void) const
{
	matrix b(row,col);
	double lambda[100];
	int counts=0,c=-1,v[100];
	ofstream fout("../inp/givens.txt");
	if(this->symmetric()!=1)
	{
		cout<<"\n Matrix is not symmetric";
		return ;
	}

	b=*this;

	for(int i=0;i<row;i++)
		for(int j=i+2;j<col;j++)
		{
			cout<<"\n Iteration   : "<<++counts;
			cout<<"\n Element : "<<b.val[i][j];
			
			if(fabs(b.val[i][j])>0.0001)
			{
				b.makezero(i,j);
				b.checktol();
			}

		}

	cout<<"\n Tridiagonal matrix is...\n";
	cout<<b;

	this->gerschgorin(lambda,c);

	b.sturm_sequence(fout,lambda,c);
	

}

